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Abstract 


The emerging field of nanomechanics is providing a new focus in the 
study of the mechanics of materials, particularly in simulating 
fundamental atomic mechanisms involved in the initiation and evolution of 
damage. Simulating fundamental material processes using first principles 
in physics strongly motivates the formulation of computational multiscale 
methods to link macroscopic failure to the underlying atomic processes 
from which all material behavior originates. This report gives an overview 
of the state of the art of the atomistic simulation of fracture and the 
application of concurrent and sequential multiscale methods to analyze 
damage and failure mechanisms across length scales. 

1.0 Introduction 

Classical fracture mechanics is based on a continuum description of material 
domains and fracture behavior described in terms of empirical parameters (Kic, J-R 
curves, Crack Tip Opening Angle, etc.). The emerging field of nanomechanics is 
providing a new insight into fracture processes beyond that available in continuum 
mechanics by simulating and characterizing fundamental atomic mechanisms involved in 
the initiation and evolution of damage. These mechanisms occur at length scales on the 
order of 10" to 1 0" m and include those leading to the creation of traction-free surfaces 
(e.g., atomic bond breakage) and plastic defects (e.g., dislocations, twins, stacking faults). 

At length scales below 10’ 6 m, direct observations of damage processes are 
extremely difficult to obtain, thus atomistic simulations have proven to be an invaluable 
tool for understanding the fundamental processes of fracture. Either quantum mechanics 
(ab-initio, tight-binding (TB) or density-functional theory (DFT)) methods or classical 
molecular dynamics (MD) or molecular statics (MS) methods can be used to simulate 
fundamental material processes using first principles in physics and provide an ultimate 
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understanding of deformation and fracture processes at the atomistic level. These 
predictions of material behavior at nanometer length scales promise the development of 
physics-based “bottom-up” multiscale analyses that can aid in understanding the 
evolution of failure mechanisms across length scales. 

This report is intended to give an overview of the state of the art in the atomistic 
simulation of fracture and the application of concurrent and sequential multiscale 
methods to analyze damage and failure mechanisms across length scales. 

2.0 Atomistic Simulations of Fracture at Interfaces 

Atomistic simulations in material science have a continuously increasing role in 
understanding the fundamental physics-based mechanisms of material behavior. The 
rapid growth of computational power and the continuous development of more robust and 
more efficient numerical methods have resulted in a substantial improvement in the 
accuracy and the performance of simulation models. Two decades ago (ca. 1980s), 
atomistic simulations were mostly used in qualitative predictions of material behavior 
(e.g., phase transitions of Lennard-Jones metallic systems, Ising model calculations in 
ferromagnetic crystals, Random Walk models for polymer chains). Currently (ca. 2008), 
it is possible to quantitatively predict many of the properties of a specific material with a 
given structural and chemical composition. The success of the transition from qualitative 
to quantitative predictions has made possible the development of new “bottom-up” 
approaches in modeling the mechanical properties of materials starting from the basic 
atomic level interactions. In these approaches, atomistic simulations play a key role in 
providing the basic information from the underlying atomistic processes that ultimately 
govern the macroscopic material parameters. 

2.1 Historical Development 

MD simulations are invaluable in studying fracture of interfaces, such as grain 
boundaries in metals, where continuum mechanics can not adequately represent structural 
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inhomogeneity. To illustrate the historical application of MD simulations to the study of 
fracture, selected analyses presented in the literature will be discussed. 


Some of the first attempts to qualitatively model the statics and dynamics of 
fracture at an atomic level date back to 1976 when Ashurst and Hoover (Ashurst and 
Hoover, 1976) used a two-dimensional triangular lattice of 512 mass points interacting 
with linear- force Hooke’ s-law springs as a representative model for a crystal (Figure 1). 
The system was evolved in a molecular-dynamics (MD) sense, through integrating 
Newton’s equations of motion for each mass point. The forces between the mass points 
were defined as linear functions of the point displacements up to a maximum strain of 
10% at which time the force was set to zero, representing bond breaking. Even for this 
simplest crystal model, the results exhibitied a number of interesting physical 
phenomena, such as varying crack propagation velocity and widespread damage around 
the crack tip, revealing details that can not be predicted using continuum mechanics. 



Figure 1. (a) Two-dimensional triangular lattice of 512 mass points interacting with 
linear- force Hooke’ s-law springs as a representative model for a crystal; (b) Image of a 
nanocrack nucleated due to atomic bond breaking in the triangular two-dimensional 
lattice. (From Ashurst and Hoover, 1976) 
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An early attempt to qualitatively simulate grain boundary segregation through 
crack propagation by MD was reported by Ishida et al. (1984) with the simulation of 
fracture in iron. These were the first simulations to reveal physical mechanisms of plastic 
processes, such as dislocation nucleation, developing at the crack tip. 

A pioneering systematic MD analysis relating the crystallography of a large set of 
tilt and twist grain boundaries in Cu and Au to their grain boundary (GB) energy, E , 
and cleavage properties was performed by Wolf some years later (Wolf, 1990). In this 
study, a large representative set of symmetric GBs in a two-dimensional space of the two 
characteristic parameters, the twist ((p) and tilt (0) angle, were investigated and their 
cleavage energies, E cl , were estimated. Cleavage energy is defined by an extension to the 
Griffith criteria for crack growth given by 

E cl = 2y — E gb (1) 

where y is the surface energy of the GB plane. Thus, a structure-energy correlation for all 
symmetric GBs in an fee metal was derived and related to the inter-granular fracture 
properties. In a related study, Yip and Wolf (Yip and Wolf, 1989) evaluated the 
atomistic concept for simulation of GB fracture. They presented an integrated approach 
for the study of the correlation between crystallographical and chemical interfacial 
structure and the physical properties relevant to intergranular fracture. This approach 
included the use of four related techniques: lattice statics for the determination of grain 
boundary energies; lattice dynamics for the analysis of local elastic constants; Monte 
Carlo simulation for determining solute segregation in grain boundary structure; and 
molecular dynamics simulation for modeling dynamic crack propagation. The major 
implication of these early works is that the problem of GB fracture must be considered 
from the point of view of the structure and properties of the GB, which significantly 
affect the process of crack propagation. 
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More recent MD simulation studies have exploited various aspects of the inter- 
granular crack propagation process. These include studies on the decohesion strength of 
a GB and its dependence on impurities (Golubovic et ah, 1995; Grujicic et ah, 1997; 
Gumbsch, 1 999); fracture stress and strain of the GB interface as a function of the atomic 
disorder due to the presence of vacancies and interstitials (Heino et al., 1998); dislocation 
emission from the crack tip as a function of the GB crystallography and structure 
(Hoagland, 1997; Cleri et al., 1999; Farkas, 2000); and crack propagation in 
nanocrystalline metals (Farkas et al., 2002; Rudd and Belak, 2002). 

Advances in computational technology have already made simulations of crack 
behavior in a single crystal of 10 atoms (Abraham, 1997; Zhou et al., 1997; Abraham, 
2001; Bulatov et al. 1998) or even 10 9 atoms (Abraham et al., 2002) possible (Figure 2). 



Figure 2. A billion atom molecular-dynamics simulation of a crack growing in a 
single crystal of Cu. Intensive emissions of dislocations from the crack tips are 
observed. (From Abraham et al., 2002) 
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Although they remain impractical for systematic studies, or for use by most practitioners, 
these large scale atomistic simulations are able to reveal many of the complex processes 
taking place at the crack tip, including dislocation emission and crack blunting in a three- 
dimensional environment. The atomistic simulations approach length scales where 
gradient theories of continuum mechanics (on the order of 0.1 pm) become valid, that is, 
where the discrete atomic structure of the matter starts to be smoothed out. In this way, 
length scales of the atomistic simulations have begun to overlap those of continuum 
model simulations (Bulatov et al., 1998; Cleri et al., 1998), using, for example, the finite- 
element method (FEM) (Xu and Needleman, 1994) or dislocation dynamics simulation 
(Noronha and Farkas, 2002) techniques. Still, computational demands and differing 
mechanistic paradigms give rise to various issues and limitations in pure atomistic 
simulation. 


2.2 Problematic Issues in the Atomistic Simulations 

The ability of the atomistic computer simulations to predict every atom’s position 
and velocity at any moment in time provides valuable information on fundamental 
material behavior not accessible from experiments. However, the necessity of 
representing the material atom by atom introduces several major issues that must be 
addressed. First is the issue of the extremely small length-scale, typically less than 100 
nm, which may make the system behavior a strong function of its boundary conditions. 
Second is the issue of the short timescale, of the order of pico and nanoseconds. For 
example, billion atom atomistic simulations remain severely limited in time scale; even 
the most powerful current supercomputers available today (ca. 2008) can not simulate the 
system response beyond the nanosecond range. This short timescale prevents the studies 
of low-rate or low-temperature deformation processes, such as diffusion or migration of 
defects and structural changes that take place over longer time scales. Third is the issue 
of the accuracy of the available interatomic potential functions used in the atomistic 
simulations to define the interactions between individual atoms. These potential 
functions define the overall physical properties of the simulated material so their 
accuracy is crucial for the correct representation of the system behavior. 
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2.2.1 Overcoming Length Scale Limitations 


There are a number of methods suggested in the literature to overcome the length 
scale limitations. All of them rely on defining special types of boundary conditions that 
may aid in mimicking an infinite system. The most popular is the use of periodic 
boundary conditions (PBC). In a system modeled using PBCs (Allan and Tildesley, 
1989), atoms at the boundary are made to interact with the atoms on the opposite side of 
the system as if the system has been multiplied in space. In this way, the boundary atoms 
are placed in the same surrounding environment as the interior atoms. The system, while 
remaining finite in volume, effectively has no boundaries. The PBC method works very 
well when the correlation length (the length over which one event can affect another 
event) of the system is smaller than the system size. If this is not true, the correlation 
effects are extended beyond the system size and are periodically multiplied causing 
severe artifacts such as strong distortions of the elastic fields or correlated dynamic 
oscillations in the system. Processes that produce long-range stress-strain fields, 
including cracks, exhibit these artifacts, and may be poorly handled by PBCs. 

For the case of long range elastic fields, another type of boundary condition, 
called a “flexible” or “elastic” boundary condition (EBC) has been developed (Sinclair et 
al., 1978). In a system modeled using EBCs, the elastic strain field due to a simulated 
defect in the system, such as a dislocation or a crack, is calculated analytically under 
infinite boundary conditions. The resulting displacements calculated at the simulated 
system boundary are applied by fixing the boundary atoms at precalculated displaced 
positions. In this way, the boundary atoms mimic the strain field of an infinite 
continuum. The method was commonly used for studying the structure of dislocation 
cores and the processes of dislocation core interactions. EBCs were also widely used for 
crack simulations. 

The drawback of the EBC method is that an imposed constraint on the boundary 
atoms affects the dynamics and the thermodynamics of the system. In addition, the 
precalculated atomic displacements are valid only for the initial configuration of the 
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system. As the system evolves in time, these boundary displacements increasingly 
deviate from the actual stress-strain state of an infinite system. An obvious solution to 
this problem would be to update the EBCs periodically. This was difficult to perform in 
early simulations as the elastic field of the evolved system often became impossible to 
solve analytically. Recently, numerical techniques have been developed to solve the 
elastic field of the evolved atomistic system and to update the displacements of the 
boundary atoms (Ohsawa and Kuramoto, 1999). These techniques have naturally 
evolved into multiscale coupled models where atomistic and continuum approaches are 
interconnected (see Section 3). Still, even with continuously updated displacements of 
the boundary atoms, the dynamics of the system remain artificially constrained. Thus, 
EBC types of methods are commonly considered “static” methods. 

An alternative to the EBC method is to apply external forces, instead of 
prescribed displacements, to the boundary atoms so that the stress at the boundary mimics 
the stress of an infinite system. The method is referred to as the “force (or traction) 
boundary condition” (FBC) method and was first introduced by DeCelis et al. (1983) to 
study fracture in iron and copper. The authors reported three major advantages of the 
FBC method versus the EBC method: (i) The FBC method allows for non-linear effects 
in the system as the boundary atoms remain unconstrained in their positions, while the 
EBC method uses linear elastic displacements, which lead to artifacts in the internal 
response; (ii) the thermodynamics of the system remains unaffected (for example, finite 
temperature simulations are readily achievable); and (iii) plastic deformation, caused for 
example by dislocation glide, is easily accommodated by the unconstrained boundaries in 
the FBC approach, while being restricted in the EBC approach. 

As in the EBC approach, externally applied forces in the FBC approach also need 
to be periodically updated to follow the dynamics of the atomistic system. A drawback 
of the FBC method is the difficulty in precisely calculating the applied external forces at 
the atomic level to closely reproduce the atomic interactions of an infinite system. 
Reaction forces calculated by a continuum mechanics model are too approximate, and 
spurious effects, such as surface tension forces due to broken atomic bonds, appear and 



make the system difficult to control and stabilize numerically. A method to eliminate 
surface tension was suggested by Cleri (Cleri, 2001) in the so-called “constant traction 
boundary conditions” molecular dynamics. This method of “constant traction boundary 
conditions” is currently used at NASA Langley Research Center (LaRC) to develop a 
novel multiscale coupling method based on the FBC approach and will be discussed in 
Section 3.2.8. 


2.2.2 Overcoming Time Scale Limitations 

The limitations on the length of system response time that can be simulated 
represent a serious obstacle in making useful predictions with atomistic MD simulations. 
To follow the atomic vibrations in an atomistic system, one has to integrate the atomic 
trajectories with a timestep of a few femtoseconds (10’ 15 s). Even with the current state- 
of-the-art gigahertz processors, the typical achievable simulation time is less than a 
nanosecond (10’ 9 s). Because the integration process is a sequential procedure, parallel 
computers cannot extend the simulation time, but can only make possible the simulation 
of larger systems for the same simulation time. 

The very short timescale requires the application of unrealistically high strain 
rates for atomistic simulations (typically in the range of 10 - 10 s’ ). The high strain 
rates require high loads to be applied (of the order of gigapascals - 10 9 Pa), which are 
often similar in magnitude to the theoretical shear strength of the material considered. In 
view of these unusually high numbers for both strain rate and magnitude of stress, it may 
seem surprising that MD simulations are successful in predicting and analyzing important 
physical processes and serving as a guiding tool for some experimental investigations. 
One of the main reasons for this success is the fact that the limitation of the short 
timescale is to a large extent compensated by the small length scale of these simulations. 
Small systems have much shorter response time, which means that the events take place 
much faster in an atomistic system than in an experiment. For example, a tensile uniaxial 

n 1 

strain rate of 10 s’ applied to a 10 nm system would mean that the boundaries of the 
system will be moving apart at a speed of only 0.1 m/s. 
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The issue of the high strain rates in MD simulations was recently addressed in 
high- temperature deformation simulations of nanocrystalline Pd with grain sizes of 
approximately 10 nm (Yamakov et ah, 2002). In spite of the extremely high strain rates 

n i 

(> 10 s' ), these simulations quantitatively validated the Coble-creep equation (Coble, 
1963) describing grain boundary diffusion creep in coarse-grained materials at strain 
rates of < 10’ 4 s' 1 . However, special care must be taken in low-temperature simulations to 
ensure that a process that might otherwise dominate the deformation behavior (that is, 
under experimental observation conditions) is not inadvertently suppressed, and hence 
overlooked, during the short time window to which MD simulations are inherently 
limited. 


There are several methods, the so-called “accelerated dynamics methods” 
(recently reviewed by Voter et al., 2002), that are being developed to permit a substantial 
increase in the time window of MD simulations. In one of the methods, called the 
“hyperdynamics” method, the interatomic potential between atoms is carefully modified 
to decrease the activation barriers for so-called “infrequent events”, such as vacancy 
migration governing lattice diffusion processes. Artificially decreasing the activation 
barriers dramatically increases the probability for occurrence of these “infrequent 
events”, thus accelerating the evolution of the system by several orders of magnitude. In 
another method, called “temperature-accelerated dynamics”, the temperature of the 
system is effectively increased while filtering out transitions that would not have 
occurred at the original temperature. This method also results in a dramatic increase in 
the dynamics of the system. In some systems, the overall acceleration may be as high as 
10 orders of magnitude for various processes, thus extending the effective simulation 
time from nanoseconds to milliseconds. 

A conceptually different approach is the “parallel replica method” (Voter, 1998) 
in which the original system is multiplied into many identical replicas and each one is 
evolved independently until an “infrequent event”, such as an atomic jump between two 
lattice sides, takes place in one of the replicas. Then, the system that has undergone the 
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change is multiplied and its replicas overwrite the reminder of the systems. The process 
is repeated until the next “infrequent event” takes place in one or more of the replicas. 
This approach is considered to be the most accurate approach of all because it does not 
require any modifications of the interatomic potential or the temperature of the 
simulation, which may introduce unknown artifacts in the predicted system behavior. 
Additionally, it is well suited for use in parallel computing environments. The trade-off 
for the high accuracy is the relatively small effective time acceleration factor, which is 
equal to the number of replicas used (usually equal to the number of available processors 
in a parallel computation, typically of the order of 10 -10 ). For additional time 
acceleration, the parallel replica method can be easily combined with the hyperdynamics 
method. 


2.2.3 Interatomic Potentials 

The interatomic potentials are used in atomistic simulations to define the force 
interactions between individual atoms in the simulated material. These potentials are 
analytic energy functions, which are simplified mathematical expressions that attempt to 
model the quantum mechanical interactions of electrons and nuclei. Their use is 
generally necessitated by the desire to model systems with sizes and time scales that 
exceed available computing resources required for quantum mechanics calculations. The 
goal of these potentials is to reproduce a variety of macroscopic physical properties of the 
simulated material, such as elastic constants, cohesive energy, melting temperature, 
vacancy formation energies, etc., by defining the elementary interactions between the 
atoms in a simple format suitable for fast computations. 

In a recent review article, Brenner (Brenner, 2000) has identified four critical 
properties that an analytic potential energy function must possess in order to be effective: 

1 . Flexibility: A potential energy function must be sufficiently flexible to 
accommodate as wide a range of fitting data of material properties as possible; 
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2. Accuracy: A potential should be able to accurately reproduce the relevant material 
properties; 

3. Transferability: A potential function should be able to describe at least 
qualitatively, if not with quantitative accuracy, structures not included in the 
fitting data base; 

4. Computational efficiency: Evaluation of the potential should be relatively simple 
and computationally efficient, so that it can be computed many times over for all 
atomistic interactions and be fast enough to allow for the simulation of 
sufficiently large systems. 

Satisfying all four criteria is a challenging task, and no general recipe exists for 
this purpose. While several standard potential functions have emerged for particular 
classes of materials, at present there is no definitive form that adequately describes all 
types of multi-atom bonding. Instead, potentials are often developed for specific 
applications with limited universality and transferability. 

There are several functional forms of interatomic potentials. Early atomistic 
simulations employed pair centrosymmetric potentials, such as the Lennard-Jones (LJ) 
potential, which are functions only of the relative distance between two atoms. As it was 
recognized that centrosymmetric pair potentials cannot reproduce the elastic anisotropy 
of a material (e.g., the anisotropy factor H = 2c 44 +c 12 -c u is always zero; Johnson, 1972), 
efforts were directed to incorporate three-body interaction terms or many-body 
environmental dependent terms in the analytical expressions (Brenner, 2000). 

Specifically for metals, the most successful form for a potential function is the 
“embedded-atom-method” (EAM) potential (Daw and Baskes, 1984), which is a 
combination of a centrosymmetric pairwise term, and a many-body term expressing a 
delocalized metallic bonding in the lattice. The EAM potential was found to accurately 
reproduce the properties of fee metals, but it gave unsatisfactory results for bee, hep and 
transition metals. In addition, the surface energy was systematically underestimated for 
all systems. Later developments were able to significantly improve the surface energy 
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value by increasing the interaction range of the potential and including a larger number of 
interacting neighbors (Mishin et al., 1999). 

However, non-fee materials remained a challenge until a modified embedded 
atom method (ME AM) was suggested (Baskes, 1992) and applied to the study of silicon 
and germanium. The MEAM potential was successful in representing both metallic and 
covalent bonds, which made its use in simulations of heterogeneous multicomponent 
systems possible (e.g., hydrogen embrittlement of aluminum or iron). A drawback of the 
MEAM potential is its relatively high computational complexity, especially in calculating 
the derivatives that are needed to calculate the interatomic forces. A simpler and 
computationally more efficient form of the MEAM potential - angular dependent (ADP) 
EAM, which still preserves the original MEAM functionality, has been suggested 
recently by Mishin (Mishin, 2005). At present, EAM and MEAM potentials have been 
developed for most pure metals and a number of binary alloys. 

3.0 Introduction to Multiscale Methods 

Modeling atomistic processes quickly becomes computationally intractable as the 
system size increases. With current computer technology, the computational demands of 
modeling suitable domain sizes (on the order of hundreds of atoms for quantum 
mechanics-based methods, and potentially billions of atoms for classical mechanics- 
based methods) and integrating the governing equations of state over sufficiently long 
time intervals, quickly reaches an upper bound for practical analyses. In contrast, 
continuum mechanics methods such as the finite element method (FEM) provide an 
economical numerical representation of material behavior at length scales in which 
continuum assumptions apply. This strongly motivates the development of analytical 
multiscale methods to link macroscopic failure to the underlying atomic processes from 
which all material behavior originates. 

Multiscale analysis is a class of systematic methodologies that have been 
developed to relate material behavior across length scales. Multiscale analyses attempt to 
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bridge length scales by providing different physics-based models that can appropriately 
represent damage mechanisms at each scale. In these approaches, models that best 
simulate the relevant physics at lower length scales are united with models at larger 
length scales through information transfer involving averaging, homogenization, or 
superposition schemes. The ultimate success of this approach is dependent on the 
accuracy of data linkage and the intrinsic fidelity of the physical models used. 

3.1 Classification of Multiscale Approaches 


Multiscale methods may be generally classified as either concurrent or sequential 
(Liu et al., 2004; Park et al., 2004). These methods are becoming valuable tools for 
linking micro- and macroscopic material behavior to atomic level processes. 


Concurrent multiscale analysis involves solving two or more strongly linked 
material models simultaneously. A general schematic of concurrently coupled atomic 
and continuum regions is depicted in Figure 3. In general, there exists a domain 
representing material at the atomic level (e.g. MD) and a domain representing material at 
a continuum level (e.g. FEM). An intermediate region involves an interface between 



Figure 3. General schematic of coupled atomistic/continuum regions. 
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these two different computational procedures that typically includes an overlap region of 
“pad” atoms/nodes in which different coupling schemes are used. 

Some approaches, however, can consist of a combination of both sequential and 
concurrent schemes, and others, such as a recently developed procedure based on 
multiscale boundary conditions, can fall far enough outside these two designations to 
constitute an independent approach (Park and Liu, 2004). Many multiscale modeling 
strategies have been explored in recent years, and the most well established state-of-the- 
art methods will be discussed. 

Sequential modeling typically involves some form of averaging of physical 
parameters that can serve as initial conditions or provide material parameters to another 
model which is analyzed separately. A desirable aspect of sequential methods is that 
length and time scales between independent material models do not have to be coupled. 
The averaging or homogenization of information across length scales that is inherent to 
sequential multiscale methods is depicted in Figure 4 where a notional coupling is shown 
across domains representing characteristic features at the sub-atomic through structural 
scales. 



Figure 4. Hierarchy of models over length scales. 
(Adapted from Oden et al., 2006) 
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3.2 Concurrent Multiscale Methods 


Concurrent methods are being continuously developed and enhanced by the 
materials science community. The common approach of the concurrent methods is to 
identify a small region of the simulated system where the representation of material is 
performed at the atomic level. This atomistic region is embedded into a larger 
surrounding region where the material is represented at the continuum level. In 
developing analytical approaches using this simulation paradigm, a primary concern has 
been the seamless coupling of forces and displacements between the different models at 
the interface between the two regions. The inherent mismatch between computational 
frameworks - atomistic (such as molecular-dynamics or molecular-statics) and 
continuum (such as finite element or finite difference) methods - and the differing 
representation of material properties can lead to simulation difficulties. 

In order to achieve a successful coupling, the usual approach is to refine the 
continuum representation (the FE mesh or the finite difference grid) down to the atomic 
scale by superposing each node to an atom at the interface region (Figure 4). In this way, 
the atomistic degrees of freedom - position and momentum of each atom - are identified 
directly with the continuum degrees of freedom - nodal displacements and their 
derivatives. Several extensive reviews of concurrent methods have been presented in the 
literature (Liu et al., 2004; Park and Liu, 2004; Oden et al., 2006). While numerous 
variations on basic procedures for coupling atomistic and continuum domains exist, only 
a handful have been well-developed and have gained a measure of widespread use. A 
brief review of these methods, representing the current state-of-the-art, follows. 

3.2.1 The Macroscopic, Atomistic, Ab initio Dynamics Method 

The MAAD (Macroscopic, Atomistic, Ab initio Dynamics) procedure was 
developed to simulate fracture by combining ab initio quantum analysis, molecular 
dynamics, and finite element continuum models (Broughton et al., 1999; Abraham et al., 
2000; Shen, 2004a). The ab initio analysis utilizes tight binding (TB) procedures to 
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predict bond breakage at the crack tip, molecular dynamics (MD) based on empirical 
force potentials to model the crack wake and surrounding atomic lattice, and a finite 
element (FE) model to simulate the far- field material. A total Hamiltonian describes the 
dynamics of the system by combining Hamiltonians of the three separate regions and 
their interfaces. In this approach, the FEM nodes correspond in a one-to-one manner 
with the interface atoms of the MD region. Figure 5 shows the basic division of the 
overall simulation into regions governed by different computational methods. 

Two significant issues have been raised regarding the MAAD approach. One is 
that the timestep used to integrate the governing equations in each of the three domains is 
equal to the smallest step required in any of them, causing a large increase in 
computational cost. Another issue is the lack of damping, which may be needed to 
remove spurious reflections at the interfaces between the three regions. 



Figure 5. Structure of MAAD coupling. (From Buehler, 2006) 


3.2.2 The Finite Element-Atomistic Method 

The FEAt (Finite Element- Atomistic) method (Kohlhoff, et al., 1991; Gumbsch 
and Beltz, 1995; Gumbsch, 1995) is a methodology similar to MAAD that links atomistic 
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representation to a continuum finite element field. Both FEAt and MAAD utilize a 
domain of “pad” atoms in the overlapping region of the MD-FEM interface in which 
individual atoms are directly linked to finite element nodes. These regions are shown in 
Figure 6. An interesting feature of FEAt is that non-local elasticity theory (Kroner, 1967) 
is used in the pad region to describe the continuum representing the finite range of 
atomistic forces and can be considered as a continuation of the lattice. These approaches 
have been successfully applied to the problem of crack propagation. 



Continuum domain 


Pad region 


Atomistic domain 


Figure 6. FEAt model of a crack tip in an fee crystal. 
(From Gumbsh and Beltz, 1995) 


3.2.3 The Coarse Grained Molecular Dynamics Method 

A generalized formulation of the conventional FEM utilizes FEM nodes 
superposed over the entire material domain to develop another computational scheme for 
atomistic-continuum coupling called Coarse Grained Molecular Dynamics (CGMD) 
(Rudd and Broughton, 1998, 2005). In this approach, a MD region is defined in which 
refinement is made such that a one-to-one atom-node linkage exists. Outside this region, 
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the finite element mesh is coarsened with individual nodes associated with many atoms. 
It is this coarse-grained (CG) finite element region that reduces the computational cost of 
representing the entire material domain. Thus, the MD region is solved using the 
integrated equations of motion for each atom, while the kinematics of the nodal degrees 
of freedom in the CG region are obtained using the equations of continuum FEM. A 
benefit of the CGMD method is that the interface between the MD and CG regions 
reduces spurious elastic wave scattering compared to other MD-FEM coupling methods. 
The MD and finite element CG regions are shown in Figure 7. 


Coarse-Grained Molecular Dynamics 



Figure 7. Depiction of CG and MD regions in the CGMD approach. 
(From Rudd and Broughton, 2005) 


3.2.4 The Quasicontinuum Method 

The Quasicontinuum (QC) method was originally formulated to provide a direct 
coupling of an atomistic region to a continuum domain (Tadmor et al., 1996; Miller et al., 
1998; Knap and Ortiz, 2001; Miller and Tadmor, 2002). The QC method is based on an 
atomic description of the material domain and uses deformation gradients and the 
Cauchy-Born rule (Bom and Huang, 1954) for homogeneous deformations to assign 
“representative atoms” or “repatoms” to describe individual atoms or enforce kinematic 
constraints on clusters of atoms that are considered as a local continuum. The basic 
assignment of atoms and possible partitioning of the QC region into local continuum 
domains and MD domains is shown in Figure 8. Continuum regions utilize FEM shape 
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functions over the domain, and material properties are obtained by a summation of the 
empirical potential function of the atoms contained in these regions. Remeshing is 
intermittently performed to update the model to resolve atomic scale detail where needed 
and to reform continuum atomic subdomains where deformation gradients are small to 
minimize computational cost. The method is commonly applied to 2-D problems and, 
because of the direct one-to-one relation between repatoms, was originally restricted to 
zero Kelvin temperature states. 


OOw 



Repatoms 


Local continuum 
region 


Figure 8. Repatoms used to define individual atoms and local 
continuum regions in the QC method. (From Knap and Ortiz, 2001) 


A finite temperature QC method has subsequently been developed (Shenoy et al., 
1999; Miller and Tadmor, 2002). The method offers a direct transition between atomistic 
and continuum fields and can effectively follow the evolution of atomistic mechanisms 
such as dislocation nucleation and crack propagation. However, special treatments are 
required to remove spurious “ghost forces” at the interface, to mitigate free-surface 
effects, and to account for finite temperature states. The QC method offers a general 
modeling technique and has been applied to simulate nanoindentation, fracture, 
dislocation motion, and interaction of grain boundaries. 
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3.2.5 The Bridging Domain Method 


Another representative concurrent coupling approach is the bridging domain 
method (Belytschko and Xiao, 2003; Xiao and Belytschko, 2004) and is based on an 
overlay approach in which MD and FEM representations are superposed in an interface 
region. This method relaxes the strict atom-node correspondence required in many other 
methods by allowing interpolation of FEM nodal displacements to be associated with 
atomic displacements in the bridging domain. The bridging domain with atom-finite 
element node overlap is shown in Figure 9. The method explicitly develops coupled 
energy Hamiltonians for the atomistic and continuum regions and enforces compatibility 
in the bridging domain using Lagrange multipliers. Dynamic behavior is simulated 
through an explicit algorithm that includes a multiple time-step scheme. This approach 
also avoids spurious wave reflection at the MD/FEM interface without introducing 
damping or filtering procedures. 


molecular model 1 continuum model 
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Figure 9. Application of bridging domain coupling method in 2-D. 
(From Xiao and Belytschko, 2004) 


3.2.6 The Coupled Atomistic/Discrete Dislocation Method 

The CADD (Coupled Atomistic/Discrete Dislocation) method (Shiari et al., 2005; 
Shilkrot et al., 2002a, 2002b; Curtin and Miller, 2003; Shilkrot et al., 2004) is specifically 
designed for problems in which dislocation formation and interaction are the physical 
mechanisms of interest and persist over long distances. The coupling at the interface is 
similar to the MAAD and FEAt methods in that it uses a pad region in which 
“handshaking” between atomic and continuum degrees of freedom occurs. CADD is 
formulated to identify the type of dislocation approaching the MD/FEM interface from 
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the atomic domain and then pass the dislocation into a surrounding discrete dislocation 
domain. Passing the dislocation is accomplished through an ad hoc addition and 
subtraction procedure. The operation involves adding analytical equations describing the 
displacement and stress fields of the dislocation to the continuum discrete dislocation 
domain and adding the forces associated with a negative image of the dislocation into the 
atomistic domain to realign the stacking fault, thus eliminating the original dislocation. 
The passed dislocations continue to interact with the MD domain through displacement 
boundary conditions at the interface. A detection band is created near the MD-FEM 
interface to determine the type of dislocation entering the continuum (eg., an edge or 
screw dislocation). This detection scheme is shown in Figure 10 where, for 2-D 
dislocations, active slip planes of a face centered cubic (fee) metal are separated by 
relative angles of tt/3. 
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Figure 10. Close-up of dislocation detection band 
near interface. (From Shilkrot et al., 2004) 

The continuum is assumed linear elastic such that a superposition can be made to 
decompose the continuum into an infinite domain that contains the long-range singular 
dislocation stress fields modeled using dislocation dynamics and a finite domain region 
that contains smooth displacement fields modeled using FEM. The decomposition of the 
coupled MD-FEM domain in CAAD is shown in Figure 1 1 with (I) depicting the infinite 
discrete dislocation region, (II) showing the bounded region defined by finite elements, 
and (III) illustrating the purely atomistic region. The infinite discrete dislocation region 
represents dislocations as superposed analytic solutions. In the figure, Q and df2 
represent volumes and surface regions, respectively, and T, f and u represent tractions, 
forces and displacements, respectively. The subscripts 0 , C, A, and I designate the 
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original combined system, the continuum region, the atomistic domain, and the interface, 
respectively and the and ‘ A ’ overbars represent the discrete dislocation field and the 
finite element field, respectively. 


This method directly addresses the representation of discrete dislocation (DD) 
plasticity in a continuum field using established DD methods (van der Geissen and 
Needleman, 1995). Because the dislocations are represented analytically, a 
computationally demanding full atomic simulation is not required, yielding a significant 
improvement in modeling efficiency. The transition between the atomic and continuum 
finite element domains utilizes a one-to-one node-atom linking that directly ties the 
interface atoms to the nodes in the continuum. An extension of the atomic region into the 
continuum is assumed in which pad atoms are superposed with continuum elements and 
are connected to nodal displacements (see Figure 3). These atoms minimize the effect of 
the free surface on the interface atoms but contribute a modeling error by introducing 
nonphysical stiffness along the interface. CADD is currently restricted to simulating 
dislocations in two dimensions. 
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Figure 11. Decomposition of coupled MD-FEM domain using 
linear superposition. (From Shilkrot et al., 2004) 
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3.2.7 The Equivalent Continuum Model 


A combined MD and Equivalent Continuum Model (ECM) method has been 
developed (Shen and Atluri, 2004b), which is similar to the Quasi-continuum methods, 
but uses the meshless local Petrov-Galerkin (MLPG) representation to link the MD and 
ECM regions. In general, meshless methods are developed to overcome some of the 
disadvantages of the finite element method, such as the need to interpolate discontinuous 
secondary variables across interelement boundaries and the need for remeshing in large 
deformation problems. The Cauchy-Born hypothesis is applied in the ECM region for 
determining the elastic properties of the continuum from the atomistic description of the 
system. The ECM and MD regions are depicted in Figure 12. 

As shown in Figure 12, in the MD region, the solid points represent atoms, while 
in the ECM region, the solid points represent atoms and the open points represent nodes 
used in the MLPG method. Thus, in the ECM region, atoms and nodes do not have to be 
coincident. The ECM method has been demonstrated in one-dimensional chain models 
and in the two-dimensional analysis of a one-atom thick planar graphite sheet. 


ECM Region MD Region 



Figure 12. Depiction of ECM and MD region in the MLPG 
approach to MD-FEM coupling. (From Shen and Atluri, 2004b) 
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3.2.8 The Embedded Statistical Coupling Method 


A recently developed approach to MD-FEM coupling has been developed based 
on a restatement of the typical boundary value problem used to define a coupled domain 
(Saether et ah, 2007 and Saether, et ah, 2008). The method uses statistical averaging of 
the atomistic MD domain to provide displacement interface boundary conditions to the 
surrounding continuum FEM region, which, in return, generates interface reaction forces 
applied as piecewise constant traction boundary conditions to the MD domain. The 
coupled MD-FEM regions are depicted in Figure 13. The two systems are 
computationally disconnected and communicate only through a continuous update of 
their boundary conditions. With the use of statistical averages of the atomistic quantities 
to couple the two computational schemes, the developed approach is referred to as an 
embedded statistical coupling method (ESCM) as opposed to a direct coupling method 
where interface atoms and FEM nodes are individually related. The structure of the 
ESCM approach is depicted in Figure 14 where the MD region near the FEM interface is 
partitioned into a series of Interface Volume Cells (IVCs) and Surface Volume Cells 
(SVCs) that are used to interface with the surrounding FEM nodes. The methodology is 


inherently applicable to three-dimensional domains, avoids discretization of the 
continuum model down to atomic scales, and permits arbitrary temperatures to be 
applied. 



Figure 13. An embedded MD region within Figure 14. Structure of the ESCM. 
an FEM domain. (Saether et al., 2007) (Saether et al., 2007) 
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3.2.9 Issues in the Development of Concurrent Methods 


Various modeling issues often arise that need to be addressed to formulate 
successful concurrent coupling methodologies. These issues naturally arise in the attempt 
to join computational domains that are modeled by distinctly different representational 
frameworks of the operant physics, namely, atomistic (MD) and continuum mechanics 
(FEM) methods. 

In MD, atoms are influenced by the nonlocal force interactions of neighboring 
atoms, while in typical FEM methods, all quantities are local to a material point. With 
the nonlocal interaction in the MD region and the local interaction in the continuum, 
linking the two can cause a loss of force reciprocity at the interface and a violation of 
momentum conservation due to induced “ghost” forces. From a theoretical standpoint of 
establishing the governing equilibrium equations, it has been observed (Curtin, 2007) that 
most concurrent coupling methods that obtain equilibrium from the variation of an energy 
functional lead to some sort of spurious force generation at the interface, while those 
methods based on summing forces to zero at the interface tend to avoid these spurious 
effects. 


Also, inelastic material processes, such as dislocations in the MD region, must be 
transformed into continuum plasticity in the FEM domain. This transformation is 
difficult to theorize and remains an open issue. Additionally, a common difficulty 
involves material mismatch as higher-order material constitutive behavior included in the 
energy potentials for atomic interactions is problematic to simulate in the constitutive 
relations used in the FEM. 

Finally, the MD-FEM interface typically generates spurious reflections of 
phonons that must be damped out. One possible solution involves implementation of a 
dynamic continuum to allow phonons to be transferred from the MD to the continuum 
region; however, the use of a dynamic continuum requires much smaller time steps than a 
quasi-static continuum, thereby incurring a large computational cost. Even with a 
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dynamic continuum, a coarse FEM mesh in the region far from the interface acts as a 
filter that naturally reflects stress waves of increasing wavelengths back to the MD 
region. 


3.3 Sequential Multiscale Methods 

Sequential - or “hierarchical” - methods are typically used for modelling 
heterogeneous materials for which different constitutive laws are required at different 
length scales. The total material deformation and stress and strain fields are typically 
decomposed into the sum of a coarse macro-component and a fine micro-component. 
The fine and coarse fields are determined through separate analyses and are summed to 
obtain the solution for the total field (Oden et al., 1999; Tadmor, et al., 2000; Clayton and 
Chung, 2006; Wagner and Liu, 2001; Hao et al., 2004). 

A notional flow of sequential coupling that utilizes cohesive zone models (CZMs) 
to carry information of microscopic failure mechanisms to predict damage progression at 
larger length scales is depicted in Figure 15 (from Glaessgen et al., 2005) . In Figure 15, 
CZMs provide the critical transition between inherently atomistic and inherently 
continuum representations. In the figure, & x represent individual displacements while A x 

represent relative displacements across a surface, i.e. = 5| op - Sj 301 . Elsewhere, the 
notation A x is used to represent combinations of relative displacements in mixed-mode 
applications such that A x = f (A\ , An , Am ). These quantities are often used 
interchangeably in the literature. Because of their importance and pervasive use in 
sequential multiscale analysis, the present discussion will focus on CZM-based methods. 

3.3.1 Cohesive Zone Models 

Cohesive zone models were originally developed to represent complicated 
nonlinear fracture processes in ductile and quasi-brittle materials (Dugdale, 1960; 
Barenblatt, 1962). CZMs were later developed to describe general adhesion and 
frictional slip along an interface (Maugis, 1992; Kem et al., 1998). The CZM approach 
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is formulated on a constitutive relationship based on applied tractions and relative 
displacements to represent separation in various fracture modes of two initially 
coincident surfaces. The relative displacements associated with the creation of a new 
fracture surface for the three fundamental fracture modes are depicted in Figure 16. 


Grain Decohesion and Local Opening (n) and 
Dislocation Motion Sliding (t) Modes 



a) Atomistic simulations 


b) Cohesive Zone Representation of 
Atomistic Results 



c) Cohesive Zones Embedded in 
Finite Element Models 




d) Idealization of Metallic Microstructures 


Figure 15. Multiscale analysis with cohesive zone models. 



Sliding 

motion 


Tearing motion 



MODE II 


MODE III 


Figure 16. Fundamental fracture modes in solids 
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In the CZM approach, the area under the CZM traction-displacement curve 
represents the work of separation to grow a crack in a particular fracture mode, i.e. the 
fracture toughness, and is given by Tvergaard and Hutchinson (1992) as 


A f 

G=\xdA (2) 

0 

where G c is the work of separation, r is the applied traction, A is a general form of the 
displacement, and A ! is the critical displacement at which complete separation has 
occurred and the tractions are zero. 

In general, the traction-displacement relationships t{A) are obtained by 
differentiation of a potential </>=<f>(A), which represents the free energy of decohesion. 
The selection of a potential function is typically based on recovering the assumed 
traction-displacement relationship, and particular forms are generally selected for 
analytical convenience. In practice, various forms have been used as shown in Table 1. 
The existence of a work potential yields the work of separation regardless of the shape of 
the function. 

CZMs have been developed using different work potentials and different 
mathematical forms, and applied to different material systems at different length scales. 
Because of the wide range of variability in CZM formulations and applications, a broad 
survey of representative work will be presented here. 

Cohesive properties along an interface have typically been approximated using 
empirical data to define the CZMs (Tvergaard and Hutchinson, 1992; Costanzo and 
Allen, 1995; Camacho and Ortiz, 1996; Klein and Gao, 1998; Zavattieri et al., 2001; 
Zavattieri and Espinosa, 2003; Turon et al., 2004). These models are frequently used in 
conjunction with the finite element method (FEM) to study fracture at macroscopic 
length scales in a wide variety of materials. 
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Wei and Anand (2004) have used a modified CZM model to study intergranular 
fracture in nanocrystalline Ni. In their FEM simulation, the CZM-based decohesion 
element approximated both reversible and irreversible inelastic sliding-separation 
deformations at the grain boundaries prior to failure. The parameterization of the model 
was performed by using available experimental data for stress-strain curves of 
nanocrystalline Ni with a large number of grains. Iesulauro (2002) has also applied the 
CZM technique to simulate fatigue crack initiation in A1 polycrystals in which a 
statistical representation of bulk material properties was input to the CZM. 

The shape of the CZM law represents the basic macro-scale behavior of the 
material near the crack (tip) under load. Various attempts have been made to determine 
the shape based on fundamental bonding characteristics in metals (Rose et al., 1983; 
Nguyen and Ortiz, 2002). The most commonly assumed forms of the traction- 
displacement law have been expressed as exponential, bilinear, and trapezoidal functions. 
A general review of various forms of CZMs is found in Chandra et al. (2002). Table 1 
shows a range of various CZM functions that have been presented in the literature and 
embedded into CZM elements. The table illustrates the basic form of the CZMs, their 
key parameters, and important features. 

Despite all of the forms of CZMs that have been proposed, a common 
mathematical form not shown in Table 1 is a bilinear constitutive relation. The bilinear 
form is often chosen because of its mathematical simplicity and because of its suitability 
for representing brittle and ductile fracture in metals (Yamakov, et al., 2006) and brittle 
fracture in polymeric and ceramic composite materials (Camacho and Ortiz, 1996). This 
form of CZM is shown in Figure 17. 

Bilinear types of cohesive zone models were used in several recent FE 
simulations of brittle fracture during multi-axial dynamic loading of ceramic 
microstructures (Zavattieri et al., 2001, 2003). 
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Table 1. Various cohesive zone models and their parameters. (From Chandra et al., 2002) 


Author Proposed model Model parameters Problem Model constants Comments 

(year) solved 


Barenblatt 
(1959, 1962) 



y=b y 


= y ^ (ductile) 
= (brittle) 


T ~ T 0 + Ti 

To is work of separation 
for brittle materials 
7) is work of plastic deformation 


Perfectly 

brittle 

materials 


The first to 
propose the 
cohesive zone 
concept 


Dugdale 

(1960) 




For small value of T/a y 



Yielding of 
thin ideal 
elastic-plastic 
steel sheets 
containing slits 


Plastic zone 
ranges from 
0.042 to 0.448 
(in.) 


Cohesive stress 
equated to yield 
stress of material 


Needleman ~ 
(1987) 

Polynomial fit Tj 

. • „ is work of separation 

■f Linear fit x p ... 

o are normalizing parameters 

/ ffmax is cohesive strength 

/ + 


- — , 

1 

f v T 7 ] 

u t 


Particle-matrix <1 = 10~ 9 to 
decohesion 10“ 8 m cohesive 
energy 1-10 J/m 2 
<W = 1000 

1400 MPa 
(j y = 350-450 MPa 


Phen omenol ogical 
model; predicts 
normal separation 


Rice and 

Wang 

(1989) 



Eo is initial Young’s modulus 
h is normalizing parameter 
o- n ax is maximum stress 
a is constant^ = 2y) 


Solute 

segregation 


Ascending part is 
equated to Eo 
considers normal 
separation and 
ignores shear 
separation 


Needleman 

(1990a) 

$ 

n , Exponential lit 

ST - Polynomial „ m . . 

— r , for T_. Linear 

— Exponential 11 

* >. fit for T t 

(j > x p is work of separation 
<5 are normalizing parameters 
tx ma , is cohesive strength 

Particle-matrix 

decohesion 


L 

A 




Predicts normal 
separation 


Needlcman 

(1990b) 



Exponential fit for normal T n 
Trigonometric fit for shear T, 

w\ A A 

\ \ / \t 


- V V 


<t> n . 4>< arc work of normal 
and shear separation 
S„. 6, arc critical displacements 
a m , % is cohesive strength 


Tvcrgaard 

(1990) 


Tvcrgaard 
and Hutch- 
inson (1992) 



fi„, d, are critical displacements 
» nu is cohesive strength 


r v is work of separation 
<\ is critical displacement 
ffmu is peak normal traction 
/interface strength <$|. 
arc factors governing shape 


Decohesion of 
interface under 
hydrostatic 
tension 
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whisker rein- 
forced metal 
matrix com- 
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Crack growth 
in elasto-plas- 
tic material, 
peeling of ad- 
hesive joints 
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Periodic shear 
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Pieriels shear 
stress due to slip 


Claims shape of 
separation law arc 
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S n = i, = Quadratic model 

1 x 10'* m 
E - 60 GPa 
(Young's mod) 

»y/E 0.005 
<>K»Joy - 5-9 


Xu and 

Needlcman 

(1993) 



is work of normal separation 
4>, is work of shear separation 
<$„. d| are critical displacements 
<Tiiui is cohesive strength 
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3.3.1. 1 Mixed-Mode Representation using CZMs 

Specific examples of mixed-mode applications of CZMs have been presented in 
the literature and are discussed here. Among these applications, Ortiz and Pandolfl 
(1999) developed a CZM approach for simulating 3-D crack propagation and considered 
displacement jumps associated with normal opening, S n , and shear opening, S s . They 
accounted for mode coupling by a simple device of introducing an effective opening 
displacement given by 

5 = Vp5;+ 8; (3) 

where the parameter p assigns different weights to the sliding and normal opening 
displacements. A simple model of cohesion is then obtained by assuming that the free 
energy potential, depends on £ through the effective opening displacement, i.e., 

* = <K8,q) (4) 
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where q is a collection of internal state variables that describe inelastic processes that 
coexist with decohesion. A potential of the form 


<t> = ec tA 
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is developed where e& 2.71828, <j c is the maximum normal traction, and S c is a 
characteristic opening displacement. Irreversibility is incorporated into the cohesive laws 
in the sense that the cohesive surfaces are assumed to unload linearly to the origin. The 
resulting exponential traction-displacement curve is shown in Figure 18(a-b). Similarly, 
a bilinear traction-displacement curve with a similar peak traction is shown in Figure 
18(c-d). 

A 3-D CZM accounting for Mode I opening, Mode II sliding, and Mode III 
tearing has been developed (Segurado and LLorca, 2004) to study decohesion in 
composite materials consisting of elastic spheres within an elasto-plastic matrix. The 
assumed form of the cohesive zone for normal opening is depicted in Figure 19. An 
exponential CZM is used to represent Mode I fracture while a linear relation is used for 
simulating tangential fracture modes. Unloading after softening is assumed to follow a 
path directly back to the origin. 

As mentioned previously, tractions are typically derived from a single elastic 
potential. In this case, the potential is given by 
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for A u n < A u c , where A u n is the normal relative displacement, and Au t \ and A ua are the 
tangential relative displacements between the crack surfaces. The term t c is the 
maximum normal stress carried by the interface undergoing purely normal separation, 
and A u c is the relative normal displacement at which all the cohesive forces vanish. The 
parameter y specifies the ratio of normal to shear stiffness of the interface, where y— 0 
indicates that the cohesive element only transfers normal stresses. The normal and 
tangential tractions at the interface can be computed by taking the partial derivatives of 
the potential with respect to the corresponding relative displacements. Other examples of 
similar mixed-mode CZM formulations can be found in Li and Chandra (2003), Chandra 
and Shet (2004), Li and Siegmund (2004) 



Figure 18. Two assumed CZM curves, (a) Exponential traction-displacement law with (b) 
loading-unloading rule, (c) Bilinear traction-displacement law with (d) loading-unloading 
rule. (From Ortiz and Pandolfi, 1999) 


Mixed-mode CZM models have also been determined through direct interpolation 
between individual single-mode CZMs (Turon et al., 2004). In this approach, a traction- 
displacement law is assumed in which normal and shear components of the traction and 
displacement are interpolated based on mode mixity. The interpolated or effective CZM 
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relates a combined traction to a combined displacement jump. The effective CZM is used 
to determine a state variable that describes the damage state used to modify the stiffness 
of the cohesive zone. An outline of the basic formulation follows. 



Figure 19. Normalized opening traction-displacement 
relation. (From Segurado and LLorca, 2004) 

The interpolation of CZM components is depicted in Figure 20. Relative 
displacement jumps between the upper and lower nodes, A = A top - corresponding to 
the maximum traction, indicate the onset of damage. Final failure of the cohesive zone is 
assumed after the relative displacement jump yields zero traction. 

T 



Figure 20. Interpolation of CZM components. 
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The interpolation of pure mode CZMs shown in Figure 20 is used to develop a 
mixed mode formulation in Turon et al. (2004). The coupled cohesive zone model is 
defined by displacement-based onset (A°) and final (A) criteria. This model uses an 
interpolation of “normal” and “shear” CZMs to obtain a single coupled CZM for mixed 
fracture modes. The Mode II and III displacements, Su and Sm, are combined by the root 
sum of the squares to make a combined “shear” displacement, 

= VCS U ) 2 + (5 m ) 2 (7) 

and the coupled displacement across the interface is defined by 


A = V(8i ) 2 +(8j 


( 8 ) 


where ( ) is the MacAuley bracket function, = V 2 (x + |v|). Mode mixity parameters 
are defined by 

P 2 


d s + (Sj 


B = 


1 + 2p 2 -2p 


( 9 ) 


where J3= 0 for pure Mode I loading, and J3= 1 for pure shear loading. The delamination 
onset criterion is defined by 



( 10 ) 


and the final criterion is derived from the B-K critical energy release rate expression 
(Benzeggagh and Kenane, 1996) given by 


G c = G lc + (G shearc - G Ic ) 


A 11 


' shear 


V w tot J 


( 11 ) 


where G tot = G\ + G s hear , G shear = G\\ + Gm, and 77 is an empirical factor used to correlate 
with experimental data. This leads to an expression for the final opening displacement 
jump given by 
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( 12 ) 


A f = 




B n 


A° 


f f f 

In the above expressions, S x , S u , and S m are the individual modal displacements 

corresponding to the final criterion and are user-specified. The 8 ° , , and l are the 

individual modal displacements corresponding to delamination onset and are computed 
from the user-specified initial stiffness values, k ° , k ° u , and k ° m , and peak tractions, t \ , 
t ° n , and t ° m , of the CZM. These parameters define the cohesive zone model used in 
decohesion finite elements for the continuum simulation of fracture at larger length 
scales. 


3.3. 1.2 CZM Models for Interfacial Sliding 

Cohesive zone models have been used to simulate shear resistance and friction 
during interfacial sliding, such as grain boundary sliding or sliding between two contact 
surfaces, in a material. While the process of sliding may look very similar to Mode II or 
Mode III fracture (as shown in Figure 16), there are several differences between fracture 
and interfacial sliding, making sliding a fundamentally different deformation mode. 

First, by definition, fracture is always related to the creation of a free surface, 
while in pure sliding (at the atomic scale) there is no free surface creation. For example, 
during GB sliding, which is a common deformation mechanism in creep deformation 
(Raj and Ashby, 1971; Yamakov et al., 2002) and governs superplasticity in metals 
(Ashby and Verrall, 1973), displacement takes place internally between atomic planes 
and does not require free surface creation. 

Second, while fracture is a localized process taking place at a crack tip, sliding is 
a delocalized process along the entire interface. The associated CZM curve for fracture 
has a finite displacement jump, and its integration (Eq. 2) always gives a finite energy of 
decohesion. In contrast, the relative tangential displacement between the two sliding 
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surfaces (or sliding distance) of an interface has no limit and can be as large as the size of 
the interface. The tractions (both normal and tangential) do not depend on the 
displacement but only on the properties of the interface and the applied load. The 
associated work of friction is proportional to the sliding distance and, unlike the energy 
of decohesion, is not a property of the interface. The interface dependent property is the 
friction coefficient, usually defined as the ratio between the tangential and normal 
tractions // = rl<j n (Zavattieri and Espinosa, 2003). 

Third, fracture creates strong stress gradients due to the stress intensity at the 
crack tip, while interfacial sliding between smooth surfaces preserves a uniform stress 
along the sliding planes. 

Recently, several CZM models for interfacial sliding have been used. For 
example, Wei and Anand (2004) modeled the strength of nanocrystalline nickel using 
CZM models based on experimental data from electrodeposited nanocrystalline Ni. In 
their model, the Mode I and Mode II cohesive zone curves were very similar, the only 
difference being that the Mode II curve was extended to twice as large of an opening 
displacement as Mode I while keeping a relatively small constant stiffness in the plastic 
region. This resulted in prediction of a certain amount of sliding with some hardening of 
the Ni grain boundaries before debonding. 

Warner et al. (2006) parameterized CZMs for grain boundary sliding in copper 
using atomistic simulation data. Two major assumptions were made in the 
parameterization: (i) the shear strength x cr i t of the CZM was assumed to be constant with 
shear displacement, and (ii) the shear and tensile strengths were assumed to be uncoupled 
because the authors were not able to extract any reliable data on their coupling. 
Parameterization of CZMs for sliding at contact surfaces that accounts for surface 
roughness and friction was performed by Zavattieri and Espinosa (2003) for the case of 
ceramics subjected to pressure-shear loading. 
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3.3.2 Nanoscopic Length-Scale Coupling Through CZMs 


When parameters that are typically used to define the CZMs are taken from bulk 
material properties (Chen et al., 1999; Iesulauro, 2002; Goyal et al., 2002; Camanho et 
al., 2003; Zhang and Paulino, 2005), they do not describe the actual physics of 
micro structural fracture. For metals, macroscale values of strength and toughness that are 
typically input to the CZM represent the aggregate responses of thousands or millions of 
grains, grain boundaries, and defects within the specimens from which they were 
obtained. Thus, these macroscale values do not represent the unique response of a 
particular interface at which a local fracture event might occur. 

If the microscale predictions are to become quantitative, consideration of the local 
nanoscale properties is required. One means of making this connection is to use the 
results of atomistic simulation to develop the constitutive relations of CZMs. In this 
approach, CZM representation of fracture begins at nanometer length scales in which 
atomistic simulation is used to predict decohesion based on fundamental damage 
mechanisms. These mechanisms include dislocation formation and interaction, 
interstitial void formation, and atomic diffusion. The development of these damage 
mechanisms progress into microscale processes such as local dislocation-based plastic 
deformation and small crack formation. Ultimately, damage progression leads to 
macroscopic failure modes such as distributed plastic yielding and the development of 
large cracks exhibiting Mode I, II, and III opening behavior. 

The connection of the CZM constitutive law with, albeit highly idealized, 
atomistic processes is a starting point for developing more realistic simulations that will 
eventually lead to accurate predictions of the failure properties of a large class of 
materials and microstructures, even when experimental data is not available. This 
methodology constitutes the basis for a sequential multiscale approach to damage 
modeling. 
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The first serious attempt to extract relevant parameters for CZM decohesion laws 
from atomistic (molecular-dynamics or molecular-static) simulations was made by Gall et 
al. (2000). Here, an atomistic model of an aluminum-silicon interface (Figure 21) was 
used to study interface debonding. In the simulation, the system was subject to a uniaxial 
tensile strain normal to the interface until debonding occurred. The internal stress was 
monitored and plotted as a function of the strain (Figure 22). The extracted atomistic 
stress-strain curve showed that after an initial elastic stretching, the interface begins to 
break (at 15% strain), and the stress rapidly decreases and then starts to oscillate around 
zero due to elastic spring-back effects of the layers after separation. Gall et al., also 
studied interface debonding for pure aluminum and pure silicon (Figure 23) and found 
that the Al-Si interface is weaker than either the Al-Al or the Si- Si interface. 


It is to be noted, as Gall and co-workers discuss in their paper (Gall et al., 2000), 
that such atomistically-derived stress-strain relations are not yet applicable for direct 
extraction of CZM traction-displacement relationships. The atomistic stress-strain plots 
from these simulations represent the local debonding of an ideal atomically flat interface 
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Figure 21. Atomic structure of Al-Si interface. (From Gall et al., 2000) 
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of dimensions between 10 and 80 nm that does not include the myriad variables that 
affect the interface. For example, the predicted debonding stress level of ~20 GPa for an 
Al-Si interface (Figure 23) is highly elevated compared to the experimental ultimate 
tensile strength of ~200 MPa for a cast Al-Si alloy where fractured and debonded Si 
particles were observed (Dighe and Gokhale, 1997; Samuel and Samuel, 1995). 



Figure 22. Atomistic derivation of a stress-strain relation during debonding of an Al-Si 
interface. The four snapshots above and below the plot show structural changes of the 
interface at different stages of debonding at the atomic level. (From Gall et al., 2000) 


Several other attempts to extract relevant parameters for CZM decohesion laws 
from molecular-dynamics or molecular-static simulations have been made by various 
groups in the last few years (Komanduri et al., 2001; Spearot et al., 2004). In addition to 
classical MD simulations, first principles quantum-mechanics based atomistic models 
(Raynolds et al., 1996) were also used to study adhesion in an NiAl-Cr interface. All of 
those models showed highly elevated debonding stress ranging from 15 to 50 GPa, which 
was two orders of magnitude higher than the experimentally observed strength of the 
corresponding materials. 
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Figure 23. Atomistic stress-strain relations for Al-Si, Al-Al, and Si-Si 
interfaces. (From Gall et al., 2000) 


There are various reasons for this discrepancy between atomistic simulations and 
experimental characterizations. In addition to various idealizations related to interface 
structure, an important factor is the selection of ensemble boundary conditions. The 
approach typically used is based on simulating the debonding of a perfect, flat interface 
under a constant tensile strain rate perpendicular to the interface. In these references 
(Komanduri et al., 2001; Spearot et al., 2004; Raynolds et al., 2006), the system size 
varied between 4 and 80 nm, and the dynamics of the atoms was severely constrained by 
the boundary conditions, which did not allow for Poisson lateral contraction and shear 
deformation. In addition, no stress intensity field was generated by introducing an initial 
crack to simulate a fracture response. As a result, plastic processes, such as dislocation 
slip, interface sliding, and interface diffusion, were strongly suppressed. Consequently, 
the simulated mechanism for interface decohesion in these works reproduced an idealized 
process of atomic adhesion (strength) rather than that of fracture at the interface. 

A recent methodology for extracting CZMs from atomistic simulations of crack 
propagation has been developed to more accurately represent a near-crack tip 
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configuration (Yamakov et al., 2006). The main goal of this approach was to extract and 
understand the contributions of different atomistic processes to an MD-based CZM 
decohesion law for intergranular fracture under local conditions of a propagating crack. 
The MD model used in this study was built to simulate a crack propagating under steady- 
state conditions through a flat high-energy grain boundary in aluminum. 

A new concept of defining cohesive zone volume elements (CZVE) as an 
atomistic equivalent of CZM elements was also introduced by Yamakov et al. (2006). In 
this concept, the CZM is a statistical representation of a large ensemble of CZVEs placed 
along the crack path during crack propagation (Figure 24). The resulting traction- 
displacement relationship is obtained as a statistical average of the behavior of a series of 
CZVEs placed along the crack path during crack propagation. When the propagating 
crack passes through the CZVEs during the simulation, the resulting stress and opening at 
each CZVE is determined, and profiles of the stress and opening displacement along the 
crack are produced (Figure 25). These profiles are taken at 1 ps intervals (time sufficient 
for a small, but detectable crack advance at the atomic level) during the entire process of 
crack propagation. The data from the entire set of many such profiles are collected and 
statistically averaged to produce a CZM traction-displacement relationship (Figure 26) 
characterizing the overall process of debonding of the interface under study. While the 
method is still under development, the results to date show the dependence of the CZM 
law on different plastic processes seen at the atomic scale, such as dislocation nucleation 
and twinning at the crack tip (Yamakov et al., 2006), heat dissipation (Yamakov et al., 
2007), and dynamic fracture processes (Yamakov et al., 2005). While the stress of 
debonding is still high (~5 GPa) compared to the experimental expectations, the analysis 
attempts to account for the local conditions of atomic bond breaking at a tip of an 
atomistically sharp crack and considers a variety of different nanoscale deformation 
processes inside the crack process zone. 
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Figure 24. MD simulation of an intergranular crack in aluminum. Cohesive 
zone volume elements (shown in the inset) are introduced along the grain 
boundary interface to extract a statistical CZM traction-displacement 
relationship. (From Yamakov et al., 2006). 



Figure 25. Stress and opening profiles extracted along the crack growing 
for 123 ps in a system prestressed at 4.25 GPa hydrostatic load. The 
corresponding snapshot of the crack is shown at the bottom. (From 
Yamakov et al., 2006). 
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Figure 26. Surface stress vs. crack opening curves a s yy (X) characterizing the 
propagation of the cleavage tip for three preloads. (From Yamakov et al., 2006). 


3.3.3 Decohesion Finite Element Formulations 

Once traction-displacement relationships are known via MD-based analysis or 
another means, they can be used to define effective decohesion finite elements that can 
then be placed along potential fracture interfaces in a finite element mesh. The resulting 
element formulation is referred to herein as a “decohesion” finite element. Finite element 
simulations can then be used to study failure evolution of material at larger length scales. 
The methodology is very generic and can be used for any interface that is subjected to 
fracture, whether the constitutive model is derived at the atomistic scale or determined at 
a larger length scale. 

Decohesion elements are formulated with two sets of coincident nodes defining 
two superposed surfaces that form a cohesive interface. Shape functions appropriate to 
the order of assumed variation over the element domain govern the interpolation of 
quantities over the superposed surfaces. These elements incorporate CZM laws to 
simulate opening failure modes and are placed between continuum finite elements to 
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predict softening and ultimate failure of the interface. In fracture studies of metallic 
materials at the micromechanical level, CZM elements can be used to predict 
transgranular fracture, if they are placed adaptively between continuum finite elements 
within grains, or to predict intergranular fracture, if they are placed along grain 
boundaries as shown in Figure 27. 




Figure 27. Embedding decohesion elements along GBs to study microstructural fracture. 

The decohesion element cohesive surfaces are typically defined as being initially 
coincident. Figure 28 shows the element surfaces under distortion. In the figure, and 
iT represent the upper and lower cohesive zone surfaces, respectively, represents the 
interpolated middle surface of the element geometry and is used to define the local 
element coordinate system. 



Figure 28. Decohesion finite element configuration. 
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Examples of decohesion elements are shown in Figures 29 and 30. Figure 29 
shows a 1-D 6-node decohesion element configuration. This line element is used to 
define cohesive interfaces between 2-D continuum elements. Although the element is 
quadratic, the Mode I opening direction is determined by computing an approximate 
linear gradient of the relative opening displacements, Av,-, over the element using only the 
end nodes. The degree of sliding in Mode II is computed by a summation of all relative 
displacements, A w/, over the element. The kinematics and integration scheme for this 
element may be found in Alfano and Crisfield (2001). 



Figure 29. Relative opening displacements for Mode I and Mode II 
fracture in a 1-D 6-node quadratic decohesion element. 


A 2-D 12-node decohesion element is shown in Figure 30. This element is used to 
define cohesive interfaces between 3-D continuum finite elements. The element is 
defined with respect to a local (x\y\ z 9 ) element coordinate system, and the kinematics 
follow a similar element formulation presented by Segurado and LLorca (2004). 

Various aspects of decohesion finite element formulations have been examined in 
the literature and are briefly summarized here. Convergence difficulties have been 
encountered and have been related to numerical problems caused by sharp comers in the 
CZM function during loading and unloading cycles from which large variations in the 
tangent stiffness matrix can occur (Gao and Bower, 2004). Another effect that can hinder 
convergence involves the type of numerical quadrature used to evaluate the element 
integrals. It has been found that Gaussian quadrature tends to link kinematic degrees of 
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freedom across the element, which affects the eigenmodes of the deformation states, and 
causes high traction gradients across the element domain. An improvement has been 
found in using Newton-Cotes methods, which act to uncouple the eigenmodes, thereby 
resulting in smoother traction profiles and improvement of convergence (Schellekens and 
de Borst, 1993). 



Figure 30. Relative opening displacements for Mode I and 
Mode II fracture in a 2-D 12-node quadratic decohesion 


Additional numerical aspects in the use of decohesion elements exist. To obtain 
accurate solutions, bounds on element size have been presented (Allen and Searcy, 2000; 
Tomar et al., 2004; Turon et al., 2007), and augmented solution schemes that include 
viscous damping to improve convergence characteristics have been advanced (Gao and 
de Borst, 1993). The initial slope of the traction-displacement law dictates the magnitude 
of the stiffness penalty constraint and, for macroscopic applications, follows the general 
rule that this parameter should be chosen large enough to enforce coincidence of nodes at 
the cohesive surface but not so large as to cause numerical ill-conditioning of the global 
stiffness matrix. Using a value that is less than optimally-large introduces a spurious 
compliance into the model that can alter results, particularly in models where decohesion 
elements are placed between many continuum elements to avoid an a priori assumption 
of where a dominant crack path will be generated. As the length scale decreases such 
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that the thickness of the cohesive zone can no longer be assumed infinitesimal, the initial 
stiffness is no longer a mathematical penalty parameter used to enforce a strict constraint 
on relative displacements but becomes increasingly reduced to represent the actual 
elasticity of the finite thickness cohesive zone. 

4.0 Summary 

This report has presented an overview of the current state of the art for 
understanding and predicting fracture in polycrystalline metals using principles from the 
emerging field of nanomechanics. Within this new realm in the study of the mechanics of 
materials, physics-based modeling of fracture involves the application of molecular 
dynamics simulation and methods of multiscale analysis. The combination of these 
computational frameworks, by linking damage mechanisms across length scales, 
promises future predictions of macroscopic material failure based on fundamental 
atomistic processes from which all material behavior originates. 

Atomistic simulation, often in the form of molecular dynamics, permits a direct 
simulation of atomistic mechanisms that cause fracture at the nanoscale and that 
ultimately dictate the overall strength and toughness of a material. These mechanisms 
include the stretching and breaking of the interatomic bonds and are often accompanied 
by atomic rearrangements that cause plastic deformation. However, for small submicron 
material domains requiring a billion atom simulation, the time scale remains severely 
limited; even the most powerful current supercomputers available today (ca. 2008) can 
not simulate the system response beyond the nanosecond range. In addition, the very 
short time scale requires the application of unrealistically high strain rates for atomistic 
simulations (typically in the range of 10 - 10 s" ). Therefore, because of their extreme 
computational expense, molecular dynamics simulation is not well suited to modelling 
large material domains. This limitation has necessitated the development of multiscale 
modeling methods to link atomistic simulation within a more computationally efficient 
framework. 
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Multiscale modeling methods that join two or more computational paradigms are 
often characterized as either concurrent or sequential. Concurrent multiscale modeling 
methods can be used to link molecular dynamics and a continuum mechanics method 
(e.g., finite element analysis) as a single simultaneously executed analysis. Often, the 
molecular dynamics simulation is embedded in a much larger finite element analysis. In 
this application, molecular dynamics is used to determine underlying physical processes 
and finite element analysis is used to provide appropriate far- field boundary conditions. 

Sequential multiscale modeling often uses cohesive zone models to provide the 
computational bridge between length scales. Cohesive zone models provide a 
numerically efficient means of representing a broad range of atomic-level mechanisms 
(i.e. dislocation formation and void formation) and configurations (i.e. grain boundaries) 
at length scales suitable for continuum mechanics approximations. Typically, the 
constitutive relationships used in cohesive zone models are approximated using 
macroscopically derived material properties. However, one new approach for the 
determination of the cohesive zone models extracts the constitutive relationship from 
molecular dynamics simulation. This methodology theoretically eliminates empirical 
descriptions of material degradation by deriving atomistically-based constitutive relations 
of material strength that can be applied at larger length scales to study the failure of 
polycrystal microstructures. 

The combined methods of molecular dynamic modeling and multiscale analysis 
techniques permit damage processes to be linked across length scales to ultimately 
predict macroscopic material behavior from first principles. As analysis methods 
improve, more realistic simulations will lead to better predictions of the failure properties 
of a large class of materials and microstructures, even when experimental data is not 
available or is difficult to obtain. 
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